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Abstract 

The paper deals with non-linear Poisson neuron network mod- 
els with bounded memory dynamics, that can include both Hebbian 
learning mechanisms and refractory periods. The state of a network 
is described by the times elapsed since its neurons fired within the 
post-synaptic transfer kernel memory span, and the current strengths 
of synaptic connections, the state spaces of our models being hier- 
archies of finite-dimensional components. We establish ergodicity of 
the stochastic processes describing the behaviour of the networks and 
prove the existence of continuously differentiable stationary distribu- 
tion densities (with respect to the Lebesgue measures of corresponding 
dimensionality) on the components of the state space and find upper 
bounds for them. For the density components, we derive a system of 
differential equations that can be solved in a few simplest cases only. 
Approaches to approximate computation of the stationary density are 
discussed. One is to reduce the dimensionality of the problem by 
modifying the network so that each neuron cannot fire if the number 
of spikes it emitted within the post-synaptic transfer kernel memory 
span reaches a given threshold. We show that the stationary distribu- 
tion of this 'truncated' network converges to that of the unrestricted 
one as the threshold increases, and that the convergence is at a super- 
exponential rate. A complementary approach uses discrete Markov 
chain approximations to the network process. We derive linear sys- 
tems for the stationary distributions of these Markov chains and prove 
that these distributions converge weakly to the stationary laws for the 
original processes. 
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1 Introduction 

Neurons are electrically excitable cells whose main function is to process and 
transmit information. They connect to each other to form neural networks that 
constitute core components of the nervous system, and so building and studying 
mathematical models of such networks is of key interest. To justify the modelling 
approach used in this paper (as described in detail in Section [2]), we will briefly 
describe the mechanism enabling neurons to communicate with each other. 

The anatomy of a neuron involves three distinct parts with different electrical 
activity functions: dendrites that form a tree and contain post-synaptic receptors 
(inputs), the cell body (soma) that integrates the input currents coming from the 
dendrites, and a long-limbed axon that terminates with pre-synaptic buttons (out- 
puts). A typical feature of the neuronal electrical activity is the propagation of 
membrane depolarisation. The membrane of a resting neuron is polarised. Brief 
high-amplitude depolarisations that propagate from the soma along the axon are 
called action potentials (or spikes) and have a characteristic shape. When a spike 
reaches an axonal termination that "connects" to a post-synaptic neuron, neuro- 
transmitters are released into the extra cellular space and excite receptors on the 
post-synaptic neuron (usually on dendrites). This generates a local variation of 
the membrane potential in that neuron, which propagates towards the soma. The 
soma can be seen as a spatio-temporal integrator of these post-synaptic potentials 
(PSP) to generate an output spike. The soma potential often remains close to the 
resting value for a few milliseconds after firing an action potential, which is referred 
to as the refractory period (through which the neuron cannot fire again). These 
basic elements of the neuronal information processing actually depend upon many 
different mechanisms at the molecular level, such as ionic concentrations, density 
of ion channels, axonal myelination, and types of neurotransmitters (for a review, 
we refer the reader to [T]). 

There exists extensive literature on mathematical modelling of both individual 
neurons and neural networks. We refer the reader interested in neurophysiological 
principles of neuron and brain operation to [H |27] and more advanced expositions 
in |314 132] of the circuitry of the brain. A detailed (but accessible and rather 
non-technical) discussion of brain networks, covering structural, functional, and 
effective connectivity and their respective dynamics, is presented in [33] (the book 
also contains extensive bibliography of the relevant research work in the complex 
network theory). A detailed overview of the computational modelling of nervous 
systems from the molecular and cellular level, including mathematical modelling 
of adaptation and learning, is given in [13]. Monograph [23] is a systematic study 
of the relationship of electrophysiology, nonlinear dynamics, and computational 
properties of neurons. One can also mention here [25\ [To] and refer to [T7] for a 
recent review of the literature in the area. 

As the duration of an action potential is relatively short (usually less than 
1 ms for sodium-based action potentials), for modelling purposes spikes are often 
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considered to be instantaneous. Hence a neuron can be modelled using a point 
process whose intensity depends on the past activity of the neuron, its incoming 
synaptic stimulation or other mechanisms. The use of point processes, such as 
Hawkes processes, for modelling the spiking activity of neurons dates back to 
papers [H QUI 02] and made it possible to study analytically the neuronal response 
to various input stimulations (for a review, see [E])- Applications of particular 
neural non-linear point process models to real data can be found in [28^ [29j [301 13^1 

In Section [2] we present descriptions of two network models we are dealing 
with in this paper. Both are "assembled" of non-linear Poisson neurons that 
can be viewed as extended versions of self-exciting Hawkes point processes, the 
difference between the two models being that the former has constant strengths 
of synaptic connections between neurons, whereas in the latter, to model Hebbian 
learning, we allow the strengths to change depending on the order in which the 
connected neurons are firing. Instead of using the formalism of point processes 
(as e.g. in [6j [5]), we choose an alternative description in terms of multivariate 
Markov processes whose states represent networks' spiking histories, with state 
spaces being products of hierarchies of simplices. This approach proves to be rather 
convenient and allows one to demonstrate ergodicity of the network processes 
under rather general conditions (Section [3]) and, moreover, to study the stationary 
distributions thereof. We show that the stationary distribution of the network has 
a smooth density with respect to a natural measure on the state space and give 
upper bounds for the components of that density on different components of the 
state space of the process. 

In Section [4] we discuss a way to reduce the dimensionality of the model and 
approximate its stationary distribution with more tractable objects. The approach 
is based on "truncating" the original process by "forbidding" neurons to fire once 
they have fired a given number n of spikes recently (within the "memory window" 
of the neuron). Thus modified process will still be Markovian and ergodic, its 
stationary distribution confined to a space of lower dimensionality and approxi- 
mating that of the original process at a super-exponential rate in n. In fact, such 
dynamics do make physical sense when the existence of the refractory period is 
taken into account, but one can further simplify the model by choosing an even 
lower threshold n. 

Section [5] deals with the problem of computing the stationary distributions 
of the networks. We derive systems of differential equations for the stationary 
distributions (unfortunately, they seem to be tractable in the simplest cases only, 
that are discussed as examples). Moreover, we prove that the stationary distri- 
bution of our network process can be approximated by those of discrete Markov 
chains constructed as discretised (in both time and space) versions of the process. 
Computing the stationary distributions for the chains is more feasible, as it only 
requires solving systems of (a large number of) linear algebraic equations. 
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2 Network dynamics and its description by 
Markov processes 

We consider a model neural network consisting of N neurons and M external 
sources. Both external sources and neurons can fire spikes, which are assumed to be 
generated by a random mechanism. External sources are assumed to fire according 
to independent Poisson processes with constant rates p k , k G {1,2, ...,M} (all 
the quantities related to external sources will be labelled with hats, and all the 
processes in the paper will be assumed to be right-continuous), whereas a neuron's 
instantaneous firing rate is determined by the value of the activation function of the 
so-called synaptic influx. For neuron i, the latter is the sum of all PSPs generated 
by spikes arriving to the ith neuron's synapses from external sources and other 
neurons, and also background activity. 

More precisely, assuming that {Tk t n}n£Z are times at which external source 
k G {1, . . . , M} fired, {T jtn } n^z are times at which neuron j G {1, . . . , N~\ fired, 
and iik(t) and ejj(i) are post-synaptic response kernel functions describing the 
effects on neuron i potential from accepting spikes through synapses connecting 
source k to neuron i and neuron j to neuron i, respectively, the total time t synaptic 
influx for neuron i is given by 

Ji(t) := Vi + ^2 W ik(Tk,m)£ik(t - T k)m ) + ^2 w ij( T j,n)eij(t - Tj >n ), 
k,m j,n 

where Vi = const represents the background activity for neuron i, the synaptic 
weights Wik{t) and Wij(t) can be positive (excitatory synapse) or negative (in- 
hibitory synapse) and, to reflect brain plasticity (e.g. to model Hebbian learning), 
they can depend on time as well. If, at time t, there is no synaptic connection of 
external source k to neuron i, we simply have Wn-(t) = 0, and likewise for network 
neurons' connections. 

The kernels £ik(t) > and eij(t) > are assumed to vanish outside a compact 
interval: for any k < M and i,j < N, 

e ik (t) = eij(t) = for t <£ [0, G], 9 = const > 0, 

which ensures causality and also means that the direct effect of any given spike on 
a neuron completely disappears within a finite time (for real life neurons, the 
order of magnitude of is 10 2 ms). 

In a previous series of papers by one of the authors ( |18l 120]). the case where 
all kernels were identical to some function e, but incorporating individual synaptic 
delays di k and d^: 

hk{t) = e(t - d ik ) and = e(t - d^) 

was considered. The delays account for both the axonal propagation of action 
potential up to the synaptic site, and for the diffusion time of the neurotransmitters 
in the synaptic cleft. In the present paper, we allow each synapse to have individual 
properties. 
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The effect of the synaptic influx on the behaviour of neuron i is expressed 
via an activation function Si(-), which is assumed to be continuous non-decreasing 
(and usually "S-shaped"), with 

< £- < Si(x) < Si < s := maxqj < oo, x G R, i < N. 

j 

Namely, denoting by Tj(i) the time of the last spike fired by neuron i prior to time 
t, by the cr-algebra generated by the evolution of our system up to time t and 
setting A(t) := (t, t + A) for A > 0, we have, as A ->■ 0, 

f neuron i fires during A(i) | & t ) = Si(Ji(t))r(t - T t (t))A + o(A), (1) 

where we used a left-continuous function r(-) G [0, 1] to model the existence of the 
so-called absolute refractory period, i.e. the time period during which a just fired 
neuron is unable to fire again. One can take e.g. 

r(s) := l(s £ (0, Sar]), $AR = const > 0, (2) 

the indicator function of the complement of the interval (0, Sar] (for real life neu- 
rons, 5ar is about 1 ms). Whatever the shape of r, we always assume that r(s) = 1 
for s > 0. 

In addition to ([I]), we assume that 

P(more than one neuron fires during A(t)\ J^t) = °(A)> 

which basically means that, given the past history j^t, the instantaneous firing of 
different neurons is driven by independent random mechanisms. 

Note that the widely studied classical Hawkes process (see [221 EJ [HJ [16]) 
corresponds to the identity activation function q in ([T]), and that the positivity of 
Si means that neurons can fire spikes in the absence of any external stimulation. 
The use of non-linear bounded functions Si is motivated by the experimentally 
observed saturation of the neuronal firing rate when its excitation increase. Ob- 
serve also that the temporal spread of the synaptic responses (modelled by e^) 
induces specific temporal correlations between the neuronal spike trains, which 
can be evaluated for the case of linear activation functions Si |22} \19\ 120], 

We will consider two types of models that differ in their assumptions concerning 
synaptic weights: 

Model I assumes that all the synaptic weights are constant: Wik(t) = Wik = 
const and Wij(t) = Wij = const for any k < M and i,j < N. 

Model II assumes a Hebbian learning mechanism in the form of spike-timing 
dependent plasticity (STDP) : if, within a short enough time interval, there are 
spikes at both pre-synaptic and post-synaptic sides of a connection, this can change 
the weight of the connection. The weight increases if the post-synaptic spike fol- 
lows the pre-synaptic one (reinforcement of the synapse), and decreases otherwise 
(depression of the synapse) . 

For simplicity we assume that, for each of the connections, the synaptic weight 
can assume finitely many values: for a common finite L, 

Wi k (t) e G ik := {g ik (l) < g ik (2) < ■ ■ ■< g ik (L)}, 
W tJ {t) G Gij := { 9ij (l) < 9ij (2) < < gij (L)}, 
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and the following discrete approximation of the STDP mechanisms discussed e.g. 
in [11]. For any i, j £ {1, . . . , N} and m € {1, . . . , L}, we have a collection of 
points 

— oo < Uij(m, m + 1) < Uij(m, m + 2) < ■ ■ ■ < Uij(m, L + 1) = 
= Uij(m, 1) < Uij(m, 2) < • • • < Uy(m, m) < oo 

In real life situations, the length <5i,w := max{iijj(m, m), \uij(m,m + 1)|} of the 
"learning window" is about 10 2 ms. We assume that <5t,iy < ©■ 

Now suppose that, for a given time t, one has Wij(t— ) = gij(m) and either 
t = or i = Tj(t) (i.e. one of the neurons i,j fired at time t). Then we put 

w ( f \ — n (js if J t = T ^ t ) and T i(*) ~ * G (u ij (m,d),u ij (m,d + 1)], 
\ t = Tj(t) and t-Ti(t) G ( Uij (m,d), Uij (m,d + 1)], 

d G {1, . . . ,L}, while otherwise the value of the weight remains unchanged. A 
similar rule applies to the weights Wi k . 

Note that the above mechanism allows one to model the emergence of new 
synaptic connections as well. Altogether, our network models provide a certain 
degree of biologically realism together with a mathematical framework that allows 
a tractable analysis. 

Having described the rules governing of our neural network, we will now present 
a Markov process model for it. Observe that, at time t, the knowledge of all the 
current synaptic weights and the times of all the spikes fired in the network within 
the time interval (t — Q,t] is all the information from the past and present that 
one needs to uniquely specify the probability distribution of the future evolution 
of the system. 

Therefore, to obtain a Markovian description of the network, we denote by 
v k {t) the number of spikes fired by external source k in the time window (t — Q,t], 
k < M. If i>k(t) = 0, then we say that source k was at the state 

X k (t) = (X kjl (t),X kj2 {t),X k>3 (t),...) = (0,0,0,...) £< 

at time t. If^(t) > 1, we set X k i(t) :=T k (t) — t + Q€ (0, 0], which is the time till 
the last spike fired by k prior to the "present" time t disappears from the moving 
window (s — @,s], s > t, and then we denote by X k ^{t) := T k (T k (t)—) — t + 6 G 
(0, 0], the time till the second last spike fired by k prior to time t disappears from 
the moving window (s — 0, s], and so on, so that in this case we always have 

X k (t) = (X k)1 (t),X k)2 (t), . . . ,X kMt) (t),0,0, . . k = l,...,M, 

with X k) x(t) > X kt2 (t) > • • • > X k i>k u\(t) > a.s. (as having two spikes at exactly 
the same time is a zero probability event). 

Likewise, the state of neuron i is described by the vector 

Xi(t) = (X iA (t),X h2 (t) . . .,X lMt) (t),0,0, ...), i = l,...,N, 

with Xi t i(t) > Xi >2 (t) > • • • > X iiU .( t \(t) > a.s., Ui(t) being the number of spikes 
fired by i during it — 0,t]. Now the complete history of spikes within the time 
window (t — @,t] is described by the vector 

Z(t) := (X(t);X(t)) := {X x {t), X 2 (t), . . .,X M (t); X x (t), X 2 (t), . . .,X N (t)). 
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For Model I, this vector will completely specify the state of the network. The 
state space for the process Z will be taken to be 



S~E M+N , where E := [j £ (n) , 

n>0 

is the union of simplices 

: = {(x 1 ,x 2 ,. ••) GK+ : 9 > si > x 2 > ■■■ > x n > 0; x n+m = 0, m > 0}, 

ra = 0, 1, 2, . . . Note that the n-dimensional simplex E^ is a face of the (n + 1)- 
dimensional one, E^ n+l \ n > 0. We will endow 5 with the product c-algebra 
y ■= <g®{M+N) ^ w h ere [ s trace of the cylindric <r-algebra on the space M N 
on E. 

For Model II, we need to specify in addition the state of the synaptic connec- 
tions. This can be done by using the matrices 

W(t) = (W ik (t))i< N , k < N , W(t) = (WijWijKN, 

in which to non-existent connections there will correspond zero entries. The new 
process Z* := (X; X; W; W) will have the state space 



s*-.=sx( n G ik ] x( n Gij) , 

AT US ]\S ' \a AS AT / 



i<N,k<AI ' i,j<N 

endowed with the natural product u-algebra that we will denote by S fi *. 

Model I dynamics. Suppose we are given an initial condition Z(0) £ S. 
Following the earlier description of the dynamics of our network, in the case of 
Model I the process Z is a piece-wise deterministic (linear) Markov process, which 
evolves for t > as follows. 

(i) Inside time intervals free of jumps, one has, for any k < M, i < N, n > 1, 

= -l(X k , n (t) > 0), ^& = -i(X i>B (t)>0). (3) 

This means that all the non-zero components of the process decay at the unit rate, 
and when the "first visible in the window" spike of, say, neuron i that occurred at 
time Ti^ n "disappears" from the moving time window (t — 0, t] at time t' = Tj ;n + 0, 
the number of positive components of Xi drops by one: ft(t') = vt{t'—) — 1. 

(ii) Given the state of the process is z = (x; x) G 5, where x = (x±, . . . , xm) 
has components x k = (x ki i, . . . , £fc, mfc , 0, 0, . . .) £ E^ mk ' with m k > 0, k < M, and 
likewise x = (x±, . . . , £Cj\r) has xi = {xi \ , ■ ■ ■ , Xi jTH > 0,0,...) e i < N, the 
instantaneous firing rate for source k is p k , and for neuron i it is given by 

Ri(z) := q ( Vi + ^2 W ik ^ fc (® ~~ £fc ' m ) 
^ k<M m>l 



+ X] lr 'J X/'' ! '0 ~ x J>)) r ( - 

j<7V n>l ' 



(4) 



7 



(iii) When source k fires (say, at time t' = Tj-.m), the only change in the state 
of Z is in the component : 

h(t') = h(t'-) + 1 

a.s. (as it is impossible to simultaneously "lose" a spike in the time window and 
acquire a new one), and the new values of the components are: 

Xk,2{t') = Xfc ; i(i'— ), 
Xk,z{t') = Xk t2 (t'—), 

Likewise, a spike fired by neuron i will mean similar changes in the component Xj. 

It is quite straightforward to write down the generator of the process Z, follow- 
ing the above description. The vector field specifying the dynamics of the process 
between jumps is piece-wise linear, and it changes its direction when, for one of 
the components X^ G E or Xi £ E, the respective integral curve running inside 
E^ n \ n > 1, hits the face _E?( n_1 ) of that simplex and then continues inside that 
lower dimensional simplex. The domain of the generator will consist of all bounded 
functions S \- > R that are path-continuous and differentiable for that vector field 
(cf. [M]). 

Model II dynamics. The trajectories of Z* will also be piece-wise deter- 
ministic (linear), with its first two components following ([3]) and the last two 
remaining unchanged between successive jumps. Jumps occur at the times when 
either external sources or neurons fire spikes, and, given the current state of the 
process is z* = (x;x;w;w) with w = (u^) and w = (wij), the instantaneous 
firing intensities are pk for source k and, instead of (j4|), 

R i( z *) '■= ft( v i + ™ik ^2 ~ ^ fc > m ) 

^ k<M m>l 

+ X] W V/~2 e ii(® ~ X j,n)j r (® ~ x i,l)- ( 5 ) 
j<N n>l ' 

for neuron i. 

When a spike is fired, the change in the components X and X is exactly 
the same as for Model I (see part (iii) of the description of its dynamics above), 
whereas the synaptic weight Wik can change when the spike was fired either by 
source k or by neuron i. As 5lw < ®> the state of Z* just prior to the spike 
completely specifies to what value the synaptic weight should change, according 
to the learning rules listed in the description of Model II. Similarly for the weights 
Wij that can change when the spike is fired by either of the neurons i and j. 
Thus we see that, in the case of Model II, Z* is a well-defined Markov continuous 
time process completely describing the dynamics of the system. Its generator will 
differ from the one for Z by the presence of terms related to jumps in the synaptic 
weights' values. 
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3 Ergodicity and the properties of stationary 
distributions 



In this section, we establish strong ergodicity of the Markov process Z* describing 
the dynamics of Model II. As Model I is a special case of the latter, this means that 
the process Z is also ergodic. The latter fact is actually an immediate consequence 
of Theorem 5 in [6] on stability of multivariate point processes with bounded 
memory dynamics (see also Theorem 6 in [6] for stability of a nonlinear multivariate 
Hawkes process with PSP transfer kernels having unbounded supports, and [26]). 
However, even in the case of Model I, our Markov process framework allows us to 
come up with much shorter and simpler proof of stability. 

Theorem 1. Under the stated assumptions for Model II, the process Z* is strongly 
ergodic: it has a unique stationary distribution ir on (S*, r y*) such that 

sup sup \P(Z*(t) € B | Z*(0) = z*) — 7r(-S) I —7- as t ^ oo. (6) 

Moreover, the convergence is exponentially fast. 

It would be most interesting to know the properties of the stationary distri- 
bution 7r. One basic fact that we can easily establish is that the distribution its 
of the first two components of (X(oo); X(oo); W (oo), W(oo)) ~ ir on (5, J?) has 
a density w.r.t. some naturally chosen measure. Moreover, we can obtain upper 
bounds for the density. 

More precisely, that natural measure on (S, S^ 1 ) is taken to be the product 
measure fi M+N , where 

fi(B) = Vn{B n ) for B= |J (B n x {0}), = (0, 0, . . .) E R+, 

n>0 n>0 

H n being the n-dimensional Lebesgue measure and B n being Borel subsets of the 
respective n-dimensional simplices 

4 n) := {(xi, ■ ■ ■ ,x n ) € Rl : G > x ± > x 2 > ■ ■ ■ > x n > 0} (7) 

that can be identified with = E^ n) x {0}. We use the convention that no is 
just the unit mass at 0. 

For (m; n) := (mi, . . . , thm', n±, . . . , njv) & 'L+ +N ', set 

k<M i<N 

similarly, E^ n ' n ^ is the product of the respective i^o-sets. 

To simplify the formulation of the next theorem, we will slightly abuse notation 
by identifying the sets £'( m;n ) with £^ m '™^ and so considering the latter as the 
components of the state space S (so that the components of fx are actually given 
on finite-dimensional spaces). 
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Theorem 2. Under the stated assumptions for Model II, if all the functions e^, 
€ij, $i and r are continuously differentiate, then, for any (m;n) G lJ^ JrN , the 
restriction of 113 to E^ n ' n ^ has a density ijj m>n w.r.t. fi admitting an upper bound 

II PT k ) (U c) exp{-eS^} < A s ™+ S ™ expj-0^} , 

k<M ' \<N ' 



where T,p := ^ k Pk, A := max{max fc p k , max; q}, S m := ^2 k < M m k , and E n 
^Zi<N n i- The density function ip m , n is continuously differentiate in the interior 
of £ , Q Tn,Tl ^ and has finite limits on its boundary. 

Remark 1. As will be easily seen from the proof of Theorem [21 if we assume 
that the function r has form ([2]) (and so is not continuously differentiable) , the 
assertion of the theorem will remain true with the only amendment that the density 
components il>m,n will be continuously differentiable inside their supports in the 
spaces of the respective dimensionalities. 

Proof of Theorem [7J It is obvious that Z* is aperiodic and stochastically continu- 
ous. So it suffices to show that the Markov "skeleton" chain Y = {Y n := Z*(Qn), 
n = 0,1,2,...}, has a recurrent state whose first hitting time distribution tail 
decays exponentially fast uniformly in the chain's initial state Z* (0) (see e.g. The- 
orem 18.1 in 0). 

First we will use a standard argument to show that the tail of r := inf{?i > 
: Z(On) = (0;0)} (i.e. the first value n such that there were no spikes in 
(0(n — l),0n]) admits such a bound. Indeed, setting for convenience P z *(-) : = 
P(-|Z*(0) = 2*), we have, for any z* G S* and t > 0, 

P z » (no spikes in (0, t]) = exp j- ^p fc + ^ • • )^dsj 

>exp{-(E^ + S ? )t}=:e-T*, (8) 

where := Yli^i, and (• • • ) represents the argument of R* along the trajectory 
of Z* on [0, t] that started at z* and experienced no jumps. 

Now setting A n := {no spikes in (0(n — l),0n]} we obtain, using recursively 
the Markov property and bound ([8]), that, for n > 0, 

P z *(t >n) ~- 




'n-l 



\m=l 



(9) 



Next we observe that (W(0n); W(@n)) is clearly an indecomposable aperiodic 
finite Markov chain, and hence it is ergodic. Take any fixed state (w';w') of this 
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chain; as it is well known, for any initial condition, the first hitting time of that 
state has an exponentially fast decaying distribution tail. Hence it is obvious that 
the state (0; 0; w ; w') £ S* will be positive recurrent for the chain Y, and that 
the tail of the first hitting time of that state by Y will admit a geometrically fast 
vanishing upper bound uniform in the initial condition of the state. The theorem 
is proved. □ 

Proof of Theorem First we will establish existence of density for transition prob- 
abilities, and then infer the desired result from that fact. 

Suppose our process started at point Z*(0) = v* and, at time 8, was at a 
point z* = (x;x;w;w) £ S* with x k £ E^ mk \ k < M, x» £ E^ Ui \ i < N. It is 
clear that the states v* and z = (x; x) £ S completely specify the trajectory of 
Z*(t) on the time interval [0, 8]; denote this trajectory by u(t), t £ [0, 8] (so that 
u(0) = v* and u(Q) = z*). Then, observing that x^i, . . . , x, jn . are the firing times 
for neuron i in the time interval [0, 8], we use the standard argument to show that 

P„* (Z(Q) £ dx\ x • • • x dxM x dx\ x ■ ■ ■ x dxjA 



(dxi) ■ ■ ■ fj, mM (dxM) 



,k<M 



e 



[] [] R t(. u (Hk-)) exp - / R*(u(t))dt 
i<N \k<m / ^ 

x fj, ni (dxi) ■ ■ {dx N ) 
--: p(v* , z)n mi (dxi) ■ ■ ■ n mu {dx M )^m (dxi) ■ ■ ■ /U njv (dx N ) 
■-p(v*,z)(i(dz). 



(10) 



Clearly, the function p(v*,z) is continuously differentiable in z in the interior of 
E^ m,n \ has finite limits on its boundary, and admits an upper bound of the form 

p(v*,z)< (J] p^Vn?r)ex P {-8E^}. (11) 
Next, in view of (|10p . for any B £ J?, we can use Fubini's theorem to write 



TT S (B) = [ 7T(dv*) [ p{v*, Z )fl(dz) 
JS* JB 



B 



7r(dv*)p(v* , z) 



n(dz) 



ip(z)fi(dz). 



(12) 



This means that tts does have density ip w.r.t. /i, and (jlip implies that %j) admits 
the desired upper bound. 

That tp is continuously differentiable in the relative interiors of the compo- 
nents of its supporting space follows from representation (|10p . the last relation 
in (|12|) and the assumption that q and the kernel functions are all continuously 
differentiable, the e's vanishing outside [0,8]. Theorem [2] is proved. □ 
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4 Approximation of tt by finite-dimensional 
distributions 



In this section we will be dealing with the simpler Model I. Even for that model, 
the state space is an infinite hierarchy of multidimensional simplices, so working 
with non-trivial distributions on it is not easy. The natural question in such a 
situation is whether one can find an appropriate approximation to the distribution 
in question, together with an approximation error bound. 

For our model, a tempting approach to finding such approximations is to con- 
sider "truncated" processes in which none of the neurons is "allowed" to 
fire more than the fixed number n > 1 times within any given time interval of 
length 0. In fact, if the model assumes existence of absolute refractory periods 
of positive length by stipulating, say, that ([2]) holds true, then that condition will 
automatically be satisfied (note, however, that one can still apply truncation with 
n < @/Sar to reduce dimensionality). 

The only difference in the dynamics of the process compared to those of 
Z is that neurons' firing intensities will now be given by 

R^iz) :=R i {z)l(x i , n = 0), i = l,...,N, (13) 

(cf. (HI). It is obvious that Z^ will also be an ergodic Markov process. Denote 
its stationary distribution on (S, 5?) by n( n \ while for the stationary distribution 
of Z (on the same measurable space) we will re- use notation ir. 

Theorem 3. Under the stated assumptions for Model I, 

sup |tt(J3) - 7r <n > (B) I < C n - (Tl+1)/2 e m , (14) 
Bay 

where C = exp{6(S /3 + S r )} and a = (1 + ln(6?))/2. 

Proof. We will use coupling. Assume that IT is a Poisson random field of unit 
intensity on IR + x R, given on some probability space, and construct a pro- 
cess {{Z{t),Z^(t))} t > with the state space S x S, whose components follow 
the original and "truncated" dynamics, respectively, start at a common state 
Z(0) = Z^(0) £ S, and are driven by the field II via the following simple mecha- 
nism. 

Introduce intervals 



m=l m=l 



k = l,...,M, 



ii:=(i?i-i,0i], i = l,...,N, 0i:=5^?j, 

i=i 

and stipulate that, in both Z and Z^ n \ external source k fires at time t if n({i} x 
Ik) > (note that {n([0, t] x Ik)}t>o are independent Poisson processes with 
constant intensities p\. , k = 1, . . . , M) . 
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Likewise, in the process Z neuron i fires at time t if 

n({t} x (^_i, + Ri{z{t-))\) > o, 

and that happens in the process Z^ at time t if 

n{{t} x (^i^-i + i??^^-))]) > o. 

Clearly, (Z, Z^) is a well-defined Markov process, and its components follow the 
desired dynamics. Note also that the process will be ergodic, like each of its 
components (this is obvious e.g. from Theorem [1]). 

Now denote by (Z(oo), Z^(oo)) a random element of S x S whose distribu- 
tion coincides with the stationary distribution of (Z,Z^). Using the standard 
argument, it is easily seen that 

\tt(B) - tt^(B)\ < P(Z(oo) ^ Z^{oo)) =: P n , Bey. 

To bound P n , denote by \Tk\, K £~N, the total length of the set 

Tr '■= {t G [0, QK] : Z(t)^Z^ n \t)} 

and observe that, from the ergodicity of (Z, Z^), one has 

Pn = Km l M- 

K^oo QK 

As we are interested in the stationary distribution, we can assume w.l.o.g. that 
the common starting point of Z and Z^ 1 ' has no components Xi in E^ m \ m > n. 
Then the trajectories Z(t) and Z^ n \t), having originated at the same point, will 
coincide with each other till the time T 1 when one of the values Xj (t) , i = 1, . . . , N, 
enters Then the respective neuron i will stay silent in Z^ at least till the 

time when the number of spikes produced by % and "visible" in the time window 
(t — 0, t] drops below n, while in Z the respective neuron will still be able to fire. 
Thus, past that time point T' , the trajectories Z(t) and Z^(t) can diverge. They 
will have to meet again, though, and the latest that will occur is at the end of the 
next "silent interval" of length 0, which, in its turn, occurs no later than at the 
time 

inf{* > T : U((t-e,t] x (-E^,E ? ] = °)}- 

To make use of the above argument to obtain an upper bound for P n , introduce 
two random sequences, {>c m }m>o an( A {7m}m>0i as follows. Letting for brevity 
9j := Qj, set 

Vij- :=n((^_i,^_i + G/2]x/i), 
V iJ+ :=n((^_ 1 + 0/2,^] x Ii), 

and then put x$ := 70 := and, for m > 1, 

x m+ i := infjj > 7 m : maxmax{T^j_, Vij+} > n/2}}, 
7 m+ i := infjj > x m+1 : n((^_i,%] x (-E^E^]) = o}. 
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Clearly, both {x m } and {7 m } are well-defined a.s. infinite increasing sequences of 
proper random variables. 

Now if Z(9 j - 1 ) = Z^(9j-i) but, for some t G {9j-i,0j] t one has Z(t) ^ 
Z( n '(t), then, for some i < N, at least one of the following two relations must 
hold: 

max{V it (j-i)-,V ii (j- 1 ) + } > -, max{y iJ _,F iJ+ } > - 

(if none of the two holds then, in any time interval of length within (Oj-2, Oj], 
neuron i will have fewer than n spikes). Thus the values >c m "mark" time intervals 
where Z and Z^ may split, whereas j m "mark" those intervals following x m 
where Z and Z^ must merge (provided that they have split indeed). 
Set 

Hk := inf{m > 1 : x m > K} - 1. 
Clearly, for t > and an arbitrary fixed e > 0, 

P(|T K | > t) < P(\T K \ > t, Hk < eK) + P(H K > eK). (15) 

First we will bound the last term. Observe that 

K 

Hk < 5^Xj> where Xj := lfmaxmax{Vjj_, Vij+] > n/2 
are i.i.d. Bernoulli random variables with success probability 



P( Xi = 1) = P(maxmax{y ij _,y iJ+ } > n/2) < 2 V P(V<j_ > n/2). 

\i<N ' * — ' 



Since has the Poisson distribution with parameter Aj := O^i/2, one can use 
Taylor's formula for the exponential series (with remainder in Lagrange form) and 
then Stirling's formula to write 

\ n / 2 i / \~( n +i)/2 

P C^-^»/2)<^f<^(5) expgd + lnA,)} 

< ^ n -(«+l)/2 e an. 



/7T 

Therefore 



p( x , = 1) < * n -(n+l)/2 e an = . 
\/7T 



Assuming that p n < 1 (otherwise the bound in the theorem will become trivial), 
and that 5 := e — p n = e — Exj > 0, we obtain 

P{H K > eK) = P(H K - KE X i > SK) 

<p(^( Xi -E Xi )>^) <e~^ K (16) 

by virtue of Theorem 10 from Chapter 5 of [3]. 
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Now we will turn to the first term on the RHS of (I15p . From the definitions 
of our random variables, it is obvious that \T%\ < O Yl m <H i^m — >tm), and so 

P(\T K \ >t,H K <eK)<v(Y, (7m " *m) > |) == Q- 

From the strong Markov property it follows that r\ m := ^m — x-m are i.i.d. geometric 
random variables, with P(r?i = k) = q(l — <7) fc_1 , A: = 1, 2, ... , where 

ff :=p(n((t-e,i] x (-S^Sf]) =0) = exp{-6(£ /; + ^)}. 

Clearly, Er/i = 1/q and <^(a) := Ee a,?1 < oo for a < — ln(l — g), 

¥>(a) = l + ^ + o(l), a^O. (17) 

Therefore, assuming w.l.o.g. that eK is integer, we have by the exponential Cheby- 
shev's inequality that 

Q < {v(a)) eK e-^ = exp{-^(^ - }. 

One can see from (]17p that, choosing t = tx •'= 0£-K"(l + /i)/<7 for an arbitrary 
fixed /i > 0, we will have, for small enough a, the bound 

Q < e~ for some c = c(a, /i). 

From here, f|15j) and (|16p we obtain the bound 

P(|Tk| >t K ) < e ~ ecK + e~ 262K . 

Clearly, ^ K P(|Tx| > iftr) < oo, and so, by Borel-Cantelli lemma, with probability 
one we have \Tk\ < tx for all large enough K. Therefore, 

e(l + h) 

P n < hmsup — — = . 

K^oo KB q 

As this holds for any e > p n and h > 0, we conclude that P n < p n /q, which 
completes the proof of the theorem. □ 



5 Computing the stationary distribution 

It is not difficult to derive differential equations (and boundary conditions for 
them) that the components ip m ,n of the stationary density of Z will satisfy in the 
case of Model I. They may be derived from the general relation 

EA/(Z(oo)) = 0, (18) 

where A is the infinitesimal generator of the process, / a function from a suitable 
subset of the domain of A, and, as before, Z(oo) ~ n. It may be easier, however, 
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to obtain them via a direct argument, making use of our Theorem [2] (of which 
the conditions will be assumed satisfied in this section unless we explicitly state 
otherwise). 

To show how to do that, we will first consider the simple case of a network 
with one external source and one neuron (with feedback). Suppose that the neu- 
ron has an absolute refractory period (so that the state space is actually finite- 
dimensional) . For simplicity, we assume throughout this section that = 1 (which 
clearly does not restrict generality). 

In this case, the state space of the process is just E x E, so that each state 
(x;x) = (xi, X2, ■ ■ ■ , x m , 0, 0, . . . ; x±, X2, ■ ■ ■ , x n , 0, 0, . . .) (note that here we sup- 
press the unnecessary first subscript indicating the number of the external source 

or neuron; likewise, p will denote here p\ etc.) belongs to one of the components 

E (m) x E {n) 

, ttx, n ^ 0. The respective density components we will denote by i\irn n* 
The first of them, V'o.O; is just the stationary probability of the silent state, for 
which we have, for 5 \ 0, 

ifa fi = p(z(S) = (o,o)) 

= P(Z(<5) = (0, 0) | Z(0) = (0, 0))P(Z(0) = (0, 0)) 

+ 



/ P(Z(6) = (0, 0) I Z(0) = (y, 0))ipi,o(v)dy 
Jo 

r 5 

+ / P(Z(<5) = (0, 0) | Z(0) = (0, y))My)dy + 0(8 2 
Jo 



= e-^^Vo.o + [\l + o(l))tl> lt0 (y)dy + f\l + o(l))^ x {y)dy + o(5), 
Jo Jo 

where the term 0(5 2 ) corresponds to the possibility that Z(0) G £'( m ) x E^> with 
m + n > 1. From the above representation we obtain that 

(p + ?(«))^o, = <Ai,o(0) + ^o,i (0), (19) 

where, using Theorem^ we put ^1,0(0) := ^i,o(0+), Vo,i(0) := ^0,1(0+). 

In the case where mn > 0, we fix a point z = (x; x) in the interior of Eq X E^ 
(see ((7J|) and set I z (5) := I&(8) x I X (,S), where 

I A (S) := (£i,fi + 5) x • • • x (x m ,x m + 5), 
Ix(5) ■= (xi,xi + 5) x • • • x (x n ,x n + 5) 

and 5 > is small enough so that I Z (S) C E^ x Eq . Using notation z + 9 for 
shifting all the components of the vector z by the same amount 6 G K and, as we 
did it before, slightly abusing notation by identifying £K m ) x E^ with E^ x E^ 1 ' , 
we have 

P(Z(<5) G / z (<5)) = P(Z(<5) G I Z(0) G / z (5) + 5) P(Z(0) G / z (<5) + 5) 

+ P{Z(8)eI z (S), Z(0) G [(/*(*) + <*) x (0,^)] x (7,(5) + 5)) 
+ P(Z(<f) G I Z (S), Z(0) G (J 4 (5) + 5) x [(4(5) + 5) x (0, <5)]) 
+ 0(5 m+n+2 ), 5\0, 
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where the last term corresponds to the possibility of Z(0) being in a space of di- 
mensionality higher than n + m + 1. Expressing the probabilities above as integrals 
of the respective density components and using Theorem we obtain the relation 




+ (1 + o(l))5 m+n+1 hV+i iri ((&, 0; y)) + ^, n +i((y; y, 0))] + 0{$ 



■m+n+2 



Subtracting from both sides the integral of ipm,n(y + 5) over I z one can then easily 
verify that the relation implies that the following differential equation must be 
satisfied: for 9 G (0, 1 — max{xi, x%}), 



Of course, the equation will hold along the whole interval formed by the intersection 
of E^ x E^ with the straight line passing through the point z and having the 
directional vector e m+n := (1, . . . , 1) G ]R m+n , the boundary condition at its right 
point being specified by the rates of transition to E^ x E^ from the state space 
components of lower dimensionalities. For example, if x\ < x\, then the right end 
point for the interval of validity of (I20p corresponds to the point where the ray 
x + 9e n , 9 > 0, hits the "right" face of E^ (the point x + 9e rn still being in 
the interior of E^). At that location, the system can only enter the component 
E^ x E^ by a jump from 

(x + 1 — xi; x* + 1 — x\) G E^ x E^ l \ where x* := (x2, x$, . . . , x n ), 

caused by a new spike fired by the neuron. Using a probabilistic argument similar 
to the one above, it is easy to see that the following must hold: 

i>m,n{z + 1 — = R(x + 1- xi;x* + 1- xi)ip mt n-i(x + l- X\]X* + 1- xi). (21) 

A similar equation will hold in the case where x\ > x\, but then the coefficient of 
i\>m-\,n on the right hand side of the respective relation will simply be p. The case 
where only one of m, n is zero is treated similarly. 

Solving equations of the form (|20p with boundary conditions (|2ip . comple- 
mented by (fT9j) and the condition that ^2 mn J ipm,nd(fi m ® p- n ) = 1, is hardly 
possible except for the simplest cases. One such case is considered in the following 
example. 

Example 1. Consider the case of a single neuron with feedback and no external 
sources. Moreover, assume that the firing rate function has the property 



gjj ?p m ,n( z + d ) = (f> + R(* + 9))lpm,n( Z + °) 

((£ + 9, 0; x + 9)) - Wn+i((* + 0; x + 9, 0)). (20) 



R{x) = for all x G E {n \ n > 2, 



(22) 
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so that there cannot be more than two spikes in any given time interval of length 
= 1 (say, due to the length of the absolutely refractory period exceeding 1/2). 
Thus the state space of the system is just E^ x E^ x E^ (which we again 
can and will identify with E^ x E^ x Eq), the density components being ip n , 
n = 0, 1, 2 (for n > 2, all i\) n = 0). 

Using an obvious notational convention, we see that an analog of f) 191) in this 
case has the form 

12(0)^ = ^(0), (23) 
while an analog of (|20p is, in the case n= 1, 

^l = R(6)M0)-M0,O), 9 €(0,1), (24) 
with the boundary condition (an analog of (|2ip ) 

Ml) = R(0)i(>o. (25) 
When n = 2, an analog of (|20p has the following form: for any y G (0, 1), 

dMy + o,o) = R{y + 0) 9)My + M) _ My + Mj 0) = 0> q g (0j 1 _ y)j 

the right-hand side of the equation being zero due to (122p . with the boundary 
condition (again an analog of (|2ip ) 



Mhi-y) = R(i-y)Mi-y)- (26) 

The last two relations immediately imply that, for any y € (0, 1), 

My + d,o) = R(i - y)Mi -y), o e (o, i- y ). 

Therefore ^(0,0) = 12(1 - 0)^.(1 - 0), so that (JMI becomes 
#i(0) 



dO 



R{6)M0)-R{l-0)Ml-0)i 9€(0,1) 



This means that the function is symmetric about the point = 1/2, so that 
ipl(0) = ?/>i(l — 0), 9 G (0,1) (hence conditions ([23j) and ([25]) are consistent) and 
the last differential equation can be re-written as 

^l = (R(6)-R(l-9))Md), 9 €(0,1). (27) 

Setting ip(9) := exp j ffj(R(y) - 12(1 - y)) dyj , we derive from <$Sty and ([27]) that 

M0) = R(0)^(9), 9 €(0,1). 

Together with (f2l?j) this completely specifies the density function ip (computing ip 
is trivial). 
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Example 2. One can also obtain a closed form solution in the case of a single 
neuron with feedback and no absolutely refractory period, but under the special 
assumption that the neuron's PSP kernel is exponential: e(t) = e~ at l(t > 0) for 
some a > 0. Of course, this violates the basic assumption that e has a finite support 
and is smooth, but such a simplified mathematical model could still provide a useful 
approximation to more realistic ones. 

It is not hard to see that in this case the dynamics of the system can be 
described by a univariate Markov process Y(t) := ^2 n e(t — T n ) > 0. Assuming 
without loss of generality that a = 1, one can see that the process Y is driven by 
the Ornstein-Uhlenbeck type equation 

dY(t) = -Y(t)dt + dZ(t), t > 0, 

Z(t) being a pure jump process with unit jumps and instantaneous jump rate 
j(Y(t)), where j(y) := <z(v + Wy) and W is the weight of the "self-connection" of 
our neuron. The infinitesimal generator A of the process Y is clearly 

Af(x) = -xf'(x) + 1 (x)(f(x + l)-f(x)), x>0; (28) 



its domain's description can be found e.g. in Proposition 4 in [21]. It is not hard 
to see that Y is ergodic (see e.g. [3]) and so has a unique stationary distribution 
that we will again denote by n. Substituting (|2"8j) into (fT8l) (with Z(oo) replaced 
by y(oo)) yields 

poo poo / fX+1 \ 

J yf{y)K{dy) = J 7 (x) (J f{y)dy\ „(dx). (29) 
Routine calculation now leads to 

roc ry 

yf'(y)n(dy) = / f'(y)b(y)dy, b(y) := / >y(x)ir(dx), 



where we used notation x + := max{x, 0} for the positive part of x. As this equation 
holds for a large enough class of functions / (see e.g. [21]) and b is continuous and 
locally bounded, we conclude that ir has a locally bounded and continuous density 
ip on (0, oo), satisfying 

rv 

ytp{y) = / j(x)tj}(x)dx, y > 0. (30) 
J(v-i) + 

This equation can be solved recursively, on intervals J n := (n, n + 1), n ^ 0. 
Straightforward calculations show that the stationary density is given by 

ip(y) := <p n (y), y £ J n , n = 0,1,2,..., 

where the ip n are found recursively as 

<Po(y) = ^(l)exp I J dxj , y £ Jo, 

ip n (y) = $ n [<p n -i](y), y G J n , n > 1, 
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where, for / defined on J n _i with a finite f(n— ), we denote by $ n [f](y), y £ J n , 
the solution <fi of the problem 

^y) = _ lk^A f{y)i y£ j ni 0( n ) = /(n-). 

The only unknown constant is just the normalizing factor that is to be de- 

termined from Jq 00 ip{y)dy = 1. At the integer points the density ip can be defined 
by continuity. 

In the general case, one can only hope to compute approximations to the 
stationary distribution of the network. One way to do that is to discretise the state 
space and approximate the differential equations for the density components ip m ,n 
discussed at the beginning of this section with respective difference equations, and 
then to solve the latter. However, although the existence of solution to the original 
system of differential equation follows from Theorem [21 establishing its uniqueness 
and also the convergence of the solutions to the systems of approximating difference 
equations presents a challenge. We will follow an alternative approach by first 
approximating the original stochastic process with a sequence of finite ergodic 
Markov chains in discrete time, and then proving convergence of their stationary 
distributions to the desired stationary distribution of Z. In the rest of the section, 
we will be dealing with our Model I, but one can easily see that analogous results 
hold for Model II as well. The only reason why we restrict ourselves to Model I 
here is that the formulation of results for the more general model is much more 
cumbersome. 

For q G N, set h = h(q) := 1/q and denote by q E the collection of all vectors 
i = (6,6, ■••,&)£ (h%) q , such that 1 > 6 > 6 > • • • > 6 > and £n+i = 
Cn+i = • • • = 6 = for some n G {0, 1, ... , q}. In particular, the null vector and 
(1, 1 — h, 1 — 2h, . . . , h) both belong to q E. Let 

q D .— q £j 

and denote by q F the "natural embedding" q S i— >• S under which the components 
of the vector £ := . . . , £ M ; £ l5 . . . , £ N ) G q S are concatenated with infinite 
strings of zeros so that, say, £j = (6,1, ■ ■ ■ ,6,<?), with the last positive component 
being £ ijni , rii < q, becomes (6,1, • • • ,6,g, 0,0, . . .) G E { ^ H \ and by q F^ 1 (B) the 
preimage of B G 5? under the mapping q F. Finally, for £ = (6,6, • • • ,6) q^i 
let 

q Ui ■= (1, (6 - h) + , (6 - h)+, . . . , - h)+) G g E. 
Now consider a Markov chain 

q Z(s) = ( q Xi(s), . . . , q X M (s); q X x (s), . . . , q X N (s)), s = 0,1,2,..., 

in the (finite) state space q S with one-step transition probabilities specified as 

follows. Given the value q Z(s) = C = (6, • • • , £m! £1, ■ • ■ > €n) e <?^> one nas the 
following transitions for the components of the vector q Z: 

X k ( s + 1) = { With P robabilit y 1 ~ h P*» ( 31 ) 

q ' \ q U% k with probability hp^, 



20 



q Xi{s + l) 



(£i ~~ h) + with probability 1 — hRi(£), 



(32) 



qU£i with probability hRi(£), 

where the operations of subtracting a scalar and taking positive parts are under- 
stood in the component-wise sense, and all the transitions occur independently of 
each other for k < M, i < N. The transitions presented as the second options on 
the right-had sides of the above relations correspond to spike firing by the respec- 
tive sources and/or neurons in the original model, and we will keep referring to 
these events as spikes in the case of the discrete model as well. 

The next theorem provides a way for numerical calculation of the stationary 
distribution 7r of our original process Z. Endow S with the topology of component- 
wise convergence and introduce the following notation. For £ = (^1,^2, • • • G 
q E, denote by 



and 



+ &,..., £n + M,0,...,0) G q E if a < 1, 

(6 + /i,£ 3 + e« + M,o,-..,o) g q E if a = 1, 



(Zi + h,& + h,...,£ n + h,h,0,...,0) G q E if a < 1, 
+ + /i,...,&, + M,0,...,0) G g £ if 6 = 1, 



possible immediate "precursors" for the state £ of a given component of the Markov 
chain q Z, i.e. the results of "inverting" transitions in ([31 \ and (|32|) . It is not hard 
to see that the states 



v & AO ■= (v &1 • • • , v &u (t M y, v ai (^), v aN (Ijv)) e 9 5, 

where a = ...,%) G {0, 1} M and a = (a\, . . . , a^) G {0, 1}^, exhaust all 
possible precursors of the state £ := . . . , £m>£i) ■ ■ ■ >£n) ^ q$ °f our Markov 
chain, and that 



p(C\v & , a (0) ■ 



n») 1(lfe=1) (i-^; 



1(4 <i) 



lk<M 



II(^(^.«(0)) 1(6=1) (i - hRi(V &ia (C))) 1($i<1) 



i<N 



are transition probabilities from those states to £. 

Theorem 4. For any q G N, t/ie Markov chain { q Z(s)} s >o is ergodic with sta- 
tionary distribution q ir = { q ir(C), C G q S} satisfying the following system of linear 
algebraic equations: 

(&;ac)e{0,l} M+N C^qS 

Moreover, as q — >■ oo ; i/ie distributions q ir o g -F _1 converge weakly to the stationary 
distribution it of Z . 



21 



Proof. That the chain q Z is ergodic is obvious since it is finite, irreducible and 
aperiodic. The system of equations that q 7r is claimed to satisfy is just an explicit 
form of the usual matrix equation q n = q ir q P for stationary probabilities, q P being 
the transition probabilities matrix of our chain. So we only need to prove the last 
claim of the theorem. 

Recall that we used to denote a "truncated version" of the process Z, of 
which the components Xj cannot take values in spaces of dimensionality higher 
than n (see Section H|) . Here we will use the same notation for a similarly "trun- 
cated" versions where the components are likewise constrained. It is easy to see 
that the assertions of Theorem [3] remains true in this case as well (with a different 
value for C). 

Denote by q Z^ a similarly truncated version of the chain q Z and observe that 
a complete analog of Theorem [3l with the same bound as in (|14p (of which the 
right-hand side does not depend on q), will hold true for that process as well. 

Now fix an arbitrary e > and choose n so large that the right-hand side 
of (|14p is less than e. That means that the stationary distributions it and 7r( n < of 
the processes Z and Z^ n \ respectively, will be e-close in total variation, and the 
same will apply to the stationary distributions q 7r and ? 7P n ' of the processes q Z 
and qZ( n ', too, so that 



sup 

g>0 



SUp \ir{B) -7T^{B)\ + SUp \ q 7T{B) - q 7T^ {B)\ 

Bay Bc„s 



< 2e. (33) 



This observation implies that it suffices to prove the claim of Theorem 0] for 
the truncated processes Z^ and q Z^ that take values in the finite-dimensional 
space M K , K := n(M + N),. To simplify notation, we will suppress the superscript 
(n) in the next two paragraphs, so that Z will mean there Z^ etc. 

To prove convergence of the stationary distributions, first assume that Z(Q) = 
qZ(0) = £ M. K and then observe that, as q — > oo, the distributions of the 
processes { q Z([qt\)} t >o weakly converge to that of {Z(t)} t >o in the Skorokhod 
space D R ic[0, oo) (see e.g. Section 5 in Chapter 2 in [T3]; by [^J we denote the 
integral part of x). This can be seen, for instance, from Theorem 2.6 in Chapter 4 
in [T3] (in fact, the purpose of the "truncation" that we did above as the first step 
in the proof was to make the state space of the processes locally compact, which is 
one of the conditions of the theorem). Indeed, extend the domain of the transition 
operator q T of the chain q Z to all bounded measurable functions / defined on the 
state space of Z (= Z^ G R K ) by setting 

q Tf(z) :=B[f( q Z(s + l))\ q Z(s) = h[qz\], 

where [y\ denotes the vector whose components are equal to the integral parts of 
the respective components of y, and let 

T(t)f(z) := E[f(Z(u + t))\Z{u) = z] , u, t > 0, 

be the transition semigroup of Z. Then the conditions of the above-mentioned 
theorem from [14] will be met provided that we show that, for any continuous 
function / on and any t > 0, one has 

hm sup| ? T^/(z) - T(t)f(z)\ = 0. (34) 
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Because of the semigroup property, it suffices to prove that convergence for t € [0, 1] 
only, which is not hard to do using representations of the form (|10|) . 

Indeed, assume that t = 1 (recall that we assumed that = 1 here); the argu- 
ment in the case t < 1 will be similar, but we will need to integrate over subspaces 
then, which makes everything even more cumbersome. Partition the component 
£,(m,n) q £ ^ e domain of integration of T(t)f(z) into cubes of edge length h with 
vertices on the grid (/iZ) s?Tl+s " (we use here notation from Theorem [2] and ignore 
incomplete cubes, i.e. the ones that intersect the "skew" faces of E^ n ' n \ as their 
contribution to the integrals will be asymptotically negligible as q — > oo). Fixing 
one of these cubes, we observe that the probability of the arrival of the chain q Z 
starting at the point x to the "left bottom" vertex of the cube after [qt\ steps will 
be given by /j Sm+s " times a product approximating the quantity p(x,z) similar 
to p(v*,z) from (|1Q[) (recall that we are dealing with Model I here, so that we 
do not need the "extended" state variable z*). Thus q T^ qt ^ f(x) will essentially 
be an integral sum approximating the integral T(t)f(x), and as the function / is 
continuous, it is a simple technical exercise to show that (I34p holds true. 

The last step in the proof is to observe that a bound of the form ([9]) will hold 
uniformly in g,n £ N for the processes {Z^ n ' (i)}t>o and { q Z^ ([qt\)}t>o as well 
(resurrecting now the superscripts (n)). Therefore there exists a t e < oo such that 
(recall that we assumed zero initial conditions for all the processes Z^ and q Z^) 
one has 



sup sup \\P(Z< n \t £ )eB)-n^(B)\ + \P( q Z< n \qt £ )eB)- q n^(B)\ 



< e. 



Now the desired assertion follows from (I33p and the weak convergence of the dis- 
tributions of q Z( n \qt £ ) to that of Z^ n \t e ) as q — > oo. Theorem H] is proved. □ 
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